Dynamical mean-field theory for bosons 
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We discuss the recently developed bosonic dynamical mean-field (B-DMFT) framework, which 
maps a bosonic lattice model onto the selfconsistent solution of a bosonic impurity model with 
coupling to a reservoir of normal and condensed bosons. The effective impurity action is derived 
in several ways: (i) as an approximation to the kinetic energy functional of the lattice problem, 
(ii) using a cavity approach, and (iii) by using an effective medium approach based on adding a 
one-loop correction to the selfconsistently defined condensate. To solve the impurity problem, we 
use a continuous-time Monte Carlo algorithm based on a sampling of a perturbation expansion in 
the hybridization functions and the condensate wave function. As applications of the formalism we 
present finite temperature B-DMFT phase diagrams for the bosonic Hubbard model on a 3d cubic 
and 2d square lattice, the condensate order parameter as a function of chemical potential, critical 
exponents for the condensate, the approach to the weakly interacting Bose gas regime for weak 
repulsions, and the kinetic energy as a function of temperature. 

PACS numbers: 71.10.Fd 
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I. INTRODUCTION 

Particles with bosonic statistics can macroscopically 
occupy a single mode at low enough temperature, even 
in the absence of correlation. This phenomenon is known 
as Bose-Einstein condensation, and leads to a variety of 
phases in strongly correlated bosonic systems. A typical 
example is 4 He, which exhibits normal, supcrfluid, solid, 
and possibly supersolid phases. Another example are di- 
lute ultracold atomic gases, which are well described by a 
Bogoliubov Hamiltonian. Strong interaction effects can 
be induced by adding lasers producing an optical lattice. 
The resulting system is an essentially clean realization of 
the Bosc-Hubbard model, which describes the competi- 
tion between hopping and on-site repulsion. Commcn- 
surability effects in the lattice lead to a phase transition 
from a supcrfluid to a Mott insulator at integer fillings 
and strong enough interaction^ Both cold gases and 4 He 
can be controlled experimentally with great accuracy and 
are virtually free of impurities and disorder. Cold atomic 
gases have the additional flexibility of tunable interaction 
strengths, and provide the freedom of changing the mass 
by choosing different alkali atoms. 

For both 4 He and ultracold bosonic gases in an op- 
tical lattice, powerful numerical methods exist for the 
strongly-interacting regime. Path integral simulations 
based on the worm algorithm can sample up to 10, 000 
bosonic atoms at T = IK for supersolid 4 He, enabling 
the study (see Ref. H for a review) of individual defects 
such as vacancies^ dislocations^ and grain boundaries^ 
For cold gases, up to 1.5 x 10 6 atoms were studied^ and 
compared directly to experiment with excellent agree- 
ment in time-of-flight images^ Such simulations involve 
a stochastic evaluation of all connected and disconnected 
diagrams occurring in a high-temperature series expan- 
sion on a finite lattice, which is an expansion in the 
hopping or kinetic energy (over temperature) around the 



atomic limits — 

The Monte Carlo simulation of fcrmionic lattice models 
suffers from the notorious sign problem, which prevents 
the study of large systems in the most interesting param- 
eter regime. A computationally tractable approximate 
method to simulate these models is dynamical mean-field 
theory (DMFT) iilri& In these calculations, the many- 
body selfenergy is approximated by all local skeleton di- 
agrams involving local propagators only, which implies 
a selfconsistent determination of the selfenergy and the 
local propagators. Non-local contributions are neglected. 
This simplification is convenient because the approxi- 
mated selfenergy can be evaluated efficiently from an ap- 
propriately defined impurity action^ By using sign-free 
(for single-site DMFT) efficient continuous-time Monte 
Carlo solvers^ one obtains the full Green's function as 
a solution to the effective impurity action in polyno- 
mial time.— ~— The simplification of the diagrammatic 
structure 1 - 1 - allows to define DMFT for arbitrary dimen- 
sions and lattice structures. A major success of DMFT is 
the understanding of the Mott metal-insulator transition 
it has provided . 13 i 15 DMFT has been extensively used 
to study model systems and - in conjunction with band 
structure techniques - to compute material properties 
for a wide range of compounds . 15 ' 20 Several extensions 
make the approximation systematic and controlled: Clus- 
ter methods 1 ^ like the dynamical cluster approximation^! 
or the cellular DMF T 22 : 23 reintroduce momentum depen- 
dence by considering multi-site impurity clusters. Meth- 
ods like DrA 2 ^ or dual fermion o 25 ' 26 systematically con- 
sider non-local diagrams beyond DMFT. In principle, a 
diagrammatic Monte Carlo evaluation of all neglected 
contributions to the selfenergy allows to estimate its ac- 
curacy, as, e.g., recently done for the Anderson localiza- 
tion problem^ 

The formulation of an analogous dynamical mean-field 
theory for bosonic lattice models has proven difficult. On 
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the one hand, from the perturbativc diagrammatic point 
of view, the idea of retaining an infinite subclass of bare 
diagrams may seem dubious since standard Fcynman dia- 
grammatic expansions in U fail notoriously for the Bosc- 
Hubbard model due to Dyson's collapse argument: In 
the complex plane, the convergence radius is zero as all 
bosons would collapse to a single point for negative inter- 
actions, with an infinite negative energy in the thermody- 
namic limit. Hence, it seems impossible to think of mean- 
ingful diagrammatic expansions in U for bosons, and only 
a few techniques are known to deal with this problem 
such as Kleinert's variational perturbation theory^ the 
cutting off of large field contributions^ 9 - or the use of a 
sequence with appropriately chosen counter terms ren- 
dering an infinite convergence radius in the absence of 
phase transitions^ In the Baym-Kadanoff approach to 
DMFT, a subset of diagrams consisting of all contribu- 
tions from local dressed propagators are retained. These 
skeleton diagrams may have radically different conver- 
gence properties than the bare series but those properties 
are essentially unknown. In the case of absolutely con- 
vergent series the bare and skeleton series are equivalent 
and physically meaningful, but due to Dyson's collapse 
argument for the bare series there is no guarantee that 
the skeleton approach will converge. 

On the other hand, keeping knowledge of the series 
provenance and by using the effective action for a single 
site (ie., restoring the action on the basis of the series), 
may still be worthwhile. In the normal phase the usual 
DMFT formalism can be applied in the same way as it 
is usually done for fermions. By using an effective action 
on a single site (which can be solved with the method of 
choice), most of the aforementioned convergence issues 
can be sidestepped. It is in fact only the occurrence of 
a condensate (which happens in the single particle chan- 
nel) in the broken symmetry phase that poses a chall eng e 
in the development of a B-DMFT formalism (cf. Refl3lf) . 
We will define our B-DMFT theory at the one-loop level 
beyond this selfconsistently defined mean-field (conden- 
sate). As shown in App. [B] our effective action for the 
impurity problem is the most general action for an im- 
purity with a broken U(l) symmetry and a one- loop cor- 
rection^ The hybridizations are then determined self- 
consistently with the underlying lattice problem. The 
approximations involved can still be understood in the 
language of diagrams, but a full interpretation in terms 
of Baym-Kadanoff functionals remains subjective in light 
of the asymptotic series expansions. We therefore prefer 
to use alternative derivation schemes in this work that 
do not rely on an expansion in the bosonic repulsion U. 

In the limit of infinite coordination number the de- 
coupling approximation for the Bose- Hubbard model be- 
comes exact (see the Appendix of Ref. [l|). This decou- 
pling approximation is recovered in our action at the 
mean-field level, since the one-loop correction vanishes. 
The original B-DMFT paper by Byczuk and Vollhardt 33 
(as well as subsequent work presented by Hu and Tong^ 4 -) 
is based on the assumption that in this limit the kinetic 



and potential energy contributions in the broken sym- 
metry phase can still be comparable to each other, in 
apparent contrast to Ref. [J. Their work postulates a 
scaling Ansatz with different scalings for condensed and 
non-condensed bosons, leading to an effective action that 
is different from the one we shall describe. The authors of 
Ref. im performed an expansion in the inverse of the co- 
ordination number around the limit of infinite coordina- 
tion number (i.e. the static mean-field result of Ref. [l|), 
but treated the condensate in a perturbative way only 
valid on the Bethe lattice. After the publication of our 
results they corrected for this and also derived a fully 
selfconsistent (and general) version of the B-DMFT for- 
malism^ 

The virtue of extending the DMFT framework to in- 
teracting Bose systems lies in the fact that certain model 
systems otherwise not amenable to bosonic simulations, 
e.g. models with frustrated interactions, or Bose-Fermi 
mixtures can now be studied numerically. The quality 
of the approximation £(fc,w) — ► £ skel (w)[Gi oc ] is system 
dependent and needs to be established on a case-by-case 
basis. We will see that the B-DMFT approximation is 
very good for the single-component Bose- Hubbard model 
in 3d, and thus may serve as a good starting point for 
more complicated systems. 

In both the fcrmionic and bosonic versions of dy- 
namical mean-field theory, the computationally chal- 
lenging part is the solution of the impurity prob- 
lem. For fermionic systems exact diagonalization; 3 ^ 
semi-analytical resummation of diagrams) 39 ! 40 quantum 
Monte Carlo j^I and numerical renormalization group^ 
methods are in wide use for single-orbital models. In 
recent years, significant progress has been made with 
the development of diagrammatic Monte Carlo impurity 
solver techniques, based on an expansion of the partition 
function in powers of the interactio n 17 ' 18 or the impurity- 
bath hybridization , 19 ' 43 allowing access to much larger 
impurity clusters, lower temperatures, and more general 
interactions^ 

In this paper we present a detailed derivation of the B- 
DMFT equations and show how the impurity-condensate 
coupling must be chosen to obtain a consistent the- 
ory. We quantify the errors introduced by the dynamical 
mean-field approximation for a system with finite coor- 
dination number by comparing with lattice Monte Carlo 
methods, and we describe the impurity solver proposed 
in Ref. |36| in such detail that the implementation of the 
method becomes straightforward. 

The paper is organized as follows: Sec. |TT] introduces 
the Bose Hubbard model, and in Sec. IIIII we derive the 
B-DMFT formalism, which is summarized in Sec. IIVI 
In section fVl we describe the diagrammatic Monte Carlo 
impurity solver. Section fVTl discusses solvable limits while 
B-DMFT results for interacting bosons on a Bethe and 3d 
simple cubic lattice are presented in Sec. IVII1 Sec. IVIIII 
provides a summary and outlook. 
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II. THE BOSE-HUBBARD MODEL 

We consider a model of spinless bosons on a lattice, 
described by the Bosc-Hubbard Hamiltonian in standard 
lattice notation, 

H = -tJ2 bjbj + |$>(^- 1 )-^E «i, (1) 

where t denotes the hopping amplitude, U the on-site in- 
teraction and \i the chemical potential. Unless otherwise 
written, our unit will be the hopping amplitude, t = 1. 

This model has three phases: (i) a normal phase at 
high temperature, (ii) a Mott insulating phase at zero 
temperature and commensurate filling for U 3> zt, and 
(iii) a supcrfluid phase occurring at low temperature. In 
the following, we will develop an effective theory which 
goes beyond static mean-field theory by including tem- 
poral (in imaginary time) fluctuations, but which is re- 
stricted to a single site and therefore contains no momen- 
tum fluctuations. 

The local terms of the action on a single site are 

f p U 
Sint = J dTb* nt (T)(d T -fi)b int {T) + —7i int (n int -l). (2) 

The subscript 'int' denotes the internal degrees of free- 
dom. This local action correctly describes the physics 
of the system in case of zero hopping (t — 0), when a 
factorization over each lattice site is exact. Hence it al- 
ready contains the low-energy physics occurring deep in 
the Mott phase. At very high temperature the action is 
also accurate because, in terms of Feynman's path inte- 
grals, the world-lines describing the evolution of particles 
in imaginary time remain almost straight lines, i.e. few 
hopping processes are present. Our task is then to replace 
this action with an effective action that also describes the 
physics deep in the supcrfluid phase (when the interac- 
tion is weak, compared to the hopping amplitude), i.e., 
in the regime where the Bogoliubov theory of the weakly 
interacting Bose gas applies. In terms of the treatment of 
the weakly interacting Bose gas (WIBG) of Ref. IH. there 
is reason to believe that such an accurate DMFT-likc ef- 
fective action exists: To leading and sublcading order (in 
the interaction), the normal and anomalous selfenergies 
are momentum independent at zero frequency, which is 
precisely the DMFT approximation. 

Note that in Ref. an explicit small symmetry break- 
ing field was added, which introduced a gap in the spec- 
trum (and removed IR divergences in leading orders). Al- 
though this violates the requirement that the spectrum 
of a superfluid should be gapless ; 45 ' 46 it was argued in 
Ref. that the leading orders of thermodynamic observ- 
ables are found on short-range distances and provided by 
Beliaev's diagrammatic technique, while for long-range 
physics (ie, the long-range wavelength fluctuations of the 
order parameter) one has to resort to Popov's hydrody- 
namic theory. Similar considerations hold for our effec- 
tive action where the gap is not fixed but depends on the 



condensate <f>. The short-range physics can be computed, 
but the long-range fluctuations are not part of the theory. 
In particular, the explicit symmetry breaking still leads 
to a finite condensate density in two dimensions at finite 
temperature (see below, sec. IVII B[) . which greatly sim- 
plifies the theory. However, while the explicit but fixed 
symmetry breaking field of Ref. 13 made the superfluid- 
normal transition first order, we expect in our case a sec- 
ond order-transition (see below, sec. IVII A2p with non- 
universal critical exponents different from static mean- 
field theory because of the dynamical corrections in the 
hybridization terms which form a one-loop correction to 
mean-field theory^ 

III. DERIVATION OF THE B-DMFT 
EQUATIONS 

In this section we present the derivation of the B- 
DMFT equations (the action and the selfconsistency re- 
lations). For completeness we present alternative (but 
equivalent) derivations of the DMFT formalism such as a 
quantum cavity reasoning (Appendix^]) and an effective 
medium approach (Appendix [Bj) in the Appendix. Here 
we implement an expansion around the atomic limit, fol- 
lowing almost literally the lecture notes by A. Georges fiL 
and consider B-DMFT as an approximation to the kinetic 
energy functional. This derivation can to a large extent 
also be found in the supplementary material accompany- 
ing our previous Letter, Ref. [36l The atomic reference 
system is interpreted as the impurity problem. We use a 
the coupling constant integration method and introduce 
source fields (Lagrange multipliers) to constrain the con- 
densate field and the connected Green's function for the 
normal bosons to their physical values. By doing so, we 
avoid the collapse arguments associated with an expan- 
sion in U mentioned in the introduction of this paper. 

A. Expansion parameter 

We introduce a parameter a £ [0,1] such that 

ffa = |£ni(ni-l)-«*X>&- (3) 

When a = 0, the atomic limit is recovered and the par- 
tition function factorizes over all sites. When a = 1 the 
full hopping is recovered, and this is the model we are 
ultimately interested in. 

B. Source fields and constraining fields 

Constraining the normal/anomalous Green's functions 
and the condensate to specified values can be done by 
introducing conjugate source fields (Lagrange multipli- 
ers) in the action. In order to constrain the condensate 
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to *J? we introduce the source field J, and analogous for 
the connected Green's function G c with source field A. 
Throughout this document we use the Nambu notation 
in which & = (</>*,(/>), J* = (J*, J) and the individual 
components of G c and A are given by 



We can then explicitly write down the grand potential 
per site (there are N s sites) which is a functional of the 
source fields and also depends on the constraining fields, 



and 



Gc(r) = 



A(r) 



G c (r) G c (t) 
G* c (r) G c (-r) 



F(-t) 2K(t) 
2K*(t) F{t) 



(4) 



(5) 



n Q [J,*,A,G c ] = - 



-L In J V[b*,b] exp | J* dr ^ b*(-d T + M )fe 4 - H a [b*,b] 
+ I drY,{J*{r)[b i {r) - &(t)] + J(t)[6J(t) - #(r)]) 



- f dr [ dr , y j F{r-r')[8b i {r)8b*{r') + G c {r-r')] 
Jo Jo i 

- f dr [ dr'Y,K{r-r')[8b*{r)8b*{r') + G* c {r-r')} 
Jo Jo 

■ J* dr dr' K* (r - r')[Sb t (r)8h (r') +G c (r-r')]\. 



Here 8b is the non-condensed part of the operator b given 
by b = (b) + 5b. 



(G) 



C. Atomic limit : impurity model 

Let us consider the atomic limit a = 0. There the 
problem is local on every site with grand potential 

I 



ft [Jo,*,A ,G c ] = - 



-L \nj V[b*,b] exp | £ dr b*{-d T + fx)^ - - 1) 

drJ2(Jo(r)[bi(r) - Mr)} + Jo(r)[6J(r) - #(r)]) 



+ I dr ( dr' y~] F (t - T')[Sbi(r)Sbl(r') + G c (r - r')] 
Jo Jo i 

+ ( dr [ dr'J2K (r-r')l8b*(r)Sbt(r') + G* c (r-r')} 
Jo Jo i 

+ £ dr dr' Kir - r')^)^') + G c (r - t')} j . 



(7) 
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From SVLq/SJq = and 6Q.q/8Jq = we obtain 

*=<b>s imp , (8) 

and from 6VL /5F = 0, 8VL /5K = 0, 8Q, /5K* = the 
relation 



G c (r) = G(r) + 



withG(r) = -(Tb^b^O))^. The expectation values 
(. . .) Simp = Tr[Te~ s ™v . . ]/Tr[Te~ s ^] are defined with 
respect to an impurity action 



1 f /3 [ p 

-- / / dTdT'6b*(T)A (r -r')5b(r') 

2 Jo Jo 

l^i I drnir) H / rfr?i(T)[?i(T) — 1] 

Jo 2 Jo 



rfrJ (r)b(r). 



(10) 



Inverting expressions ([5]) and © yields J[<I>,G C ] and 
A [$,G C ] and thus a functional T of the condensate 
and connected impurity Green's function: 

r [*,G c ] =F imp [*,G c ] 

+ / rfr[F (r)G c (r) + K (r)G* c (r) + K* q {t)G c {t)\ 
Jo 

+ j^jyr^2[J*(r)Mr) + Jo(rW l (r)}, (11) 

where -Fi mp is the free energy of the local quantum im- 
purity model. 



where uj n = 2irn/ ft are even Matsubara frequencies and 
£k is the dispersion relation of the non-interacting bosons. 

We now arrive at the formal expression for the exact 
functional T = T a — i, 



r[*,G c ] = r [*,G c ] + £[*,G c ], 



(9) with the kinetic energy functional given by 



/C[*,G C ] = J da^[#,G c ]. 



(14) 



(15) 



Requiring stationarity (8T/54>*(t) — 0, 5T/5(f)j{r) = 0) 
determines the value of the source field conjugate to the 
condensate (assume a homogeneous condensate over the 
lattice), 



Jo = zt&. 



(16) 



Since the condensate is time-independent (and taken 
real), we drop the r dependence of Jo as well. The other 
stationarity requirement (5T/SG C = 0, 5T/SG C = 0) de- 
termines the hybridization function appearing in the im- 
purity action: 



*b(r) 



SG c (t) 



SIC 



5G c {t) 



(17) 



Note that for the case z = oo, we have identically Sb = 
following Appendix A in Refill and only static mean-field 
theory exists, which is physically clear: For a thermody- 
namic condensate to be gapless, the k = component of 
the Green function should not decay in time. The infinite 
connectivity of the lattice removes any k dependence and 
in combination with bosonic statistics one then sees that 
the decoupling approximation is exact. 



D. Full model 



E. Approximation to the kinetic energy functional 



The exact functional of the (local) Green's function 
and condensate are constructed by using the coupling 
constant integration method, starting from the atomic 
limit: 



r = r Q= i = r + / da^^-. 

Jo da 



(12) 



By using the stationarity of fl (a-dcrivatives of the La- 
grange multipliers do not contribute) we get 



dT a i 

da N 



l - J drtJ^ibUrMr)) 

l — f drrtYftXWit?) + (5bt(r)6b j (r))} 
sP Ja (ij) 

2^ Tr E e kGc(k^„)k, Gc 

n,k 



(13) 



B-DMFT can now be considered as an approximation 
to the kinetic energy functional. With the single-particle 
Green's function of the Bosc-Hubbard model in the pres- 
ence of source fields and for arbitrary coupling constants, 
we can define a selfenergy 

G"(k, iu n ) = [iu n <7?, + (/i - ae k )I + A a - £ Q [k, iw™]] -1 , 

(18) 

where 03 is the Pauli matrix with ±1 on the diagonal. 

The DMFT approximation consists in replacing the 
selfenergy S Q for arbitrary a by the impurity model self- 
energy So. Hence, 

G"(k, iu n ) I B-DMFT = 

[iu n <7 3 + (fj, - ae k )I + A a - S Q=0 [iw„]] _1 
= [A a - A + G" 1 - aekl] -1 , (19) 

where we have used that the impurity selfenergy satisfies 
the Dyson equation 

S Q= o[* w n] = G 1 — G c 1 

= iu n a 3 + fjl + Aq - G" 1 (20) 
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and Gq is the bare Green's function. Summing over 
k, and using the constraint on the local lattice Green's 
function, we obtain the following relation between G c 
and the hybridization function: 

G c (iuj n ) = [ deD{e)(( - ael)- 1 = -£>(-), (21) 
J a \as 

with £ = A Q — Ao + G^ 1 . We used the non- interacting 
density of states D(e) = S(e— et) and its Hilbert- 

transform D(z) = J deD(e)(z — el) -1 . By introducing its 
inverse, D(R(g)) = g, the relation above can be inverted 



(aR(aG c ) = C = A Q - A + G,: 1 ) and yields the hy- 
bridization function as a functional of the local Green's 
function, 

A a [iw n ;#,G c ] = -G c - 1 +A [*,G c ]+ai?(aG c ). (22) 

Inserting the above relation into (fTO)) . the lattice Green's 
function expressed as a functional of G c becomes 

G°(k, iw„)| B -DMFT = (aR(aG c ) - ae^iy 1 . (23) 

Equation (|T3"|) can now be evaluated with G"(1c)|b-dmft: 



--i-Tr^e k G^k,» W „)| B -DMFT - f drt ^[^(r)^(r)] 

= -^ Tr E / deeD(e)(R(oG e )-d)- l -^jyrt^[<l>^T)Mr)] 



(24) 



r 



An explicit expression for the B-DMFT approximation 
to JC[<&, G c ] therefore reads 



K. 



B-DMFT 



G c R(aG c ) - -I 
a 



- 3jkjf*«D*;<TW,(r)j, (25 > 



where the last term reduces to —zt<p*<p for a constant, 
homogeneous condensate. 



F. Stationarity conditions 



This equation defines the B-DMFT selfconsistency con- 
dition. 

B-DMFT maps the bosonic lattice problem to a self- 
consistent solution of an impurity model, whose action 
(expressed in terms of the full operators b, therefore 
dropping the term d T ) now takes the form 



It immediately follows from Eq. (f25[) that the sta- 
tionarity condition for the condensate is unaltered in 
the B-DMFT approximation (Jo = zt<k), while the sta- 
tionarity condition for the connected Green's function 
(5T/SG C = 0, 5T/SG C = 0) reads in the B-DMFT ap- 
proximation (use R(aG) + aGR'(aG) = d a [aR(aG)] and 
the cyclical properties of the trace) , 

A [ioj n ; *, G c ]| B -dmft = -R[G c (iw n )] + G c (iw„)" 1 
= — ioj n a 3 - fil + S iin p + G" 1 . (26) 

Here we have used again the Dyson equation for the sec- 
ond equality. Applying D{.) to both sides of Eq. (|2"6) 
gives 

G c (iuj n ) = J deD(e)(iuj n a 3 + (fi - c)I - Simp) -1 . (27) 



\ [ [ dTdrV(r)Ao(T-r')b(r') 

z Jo Jo 

—fi / drn(r) + 77/ dTn{T)[n{T) — 1] 
Jo 2 Jq 

/ r p \ r 13 

-&(zt- dT'Ao(r')J J drb(r). (28) 

Assuming K = K* , tp = (p* and dropping all subscripts 
we can write the action in our final version as 



Simp — 2 



1 r P r/3 

- / drdT y b t (r)A(T-T')b(r') 

2 Jo Jo 

f U f fi 

—fi I drn(r) H / drn(r)[n(r) — 1] 

Jo 2 J Q 

-k& / drb(r) (29) 
Jo 



Simp 



where the coefficient k is given by 



k = zt — An(iw n = 0) — Ai2(iw„ = 0). (30) 
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The solution of the impurity problem yields the con- 
densate 3? (Eq. ©), the connected Green's function G c 
(Eq. ((9])) and the selfenergy S; mp of the impurity model. 
The right hand side of Eq. (|2"T)) then defines the local 
lattice Green's function, which is identified with the im- 
purity Green's function and thus allows to obtain the 
new hybridization function for the next iteration by us- 
ing Eq. (HSJ). 



IV. DMFT PROCEDURE 

Solving the single site impurity model means comput- 
ing the local Green's function 



G(r) = -<Tb(r)bt(0)) Sl , 



(31) 



and the condensate <fr for the impurity action given by 
Eq. The hybridization function A(r) and the con- 

densate order parameter which is constant in time, 
are calculated by a selfconsistent procedure starting from 
some arbitrary initial values (usually obtained from the 
non-interacting or static mean- field limit). We then solve 
the impurity problem and calculate the new updated pa- 
rameters A(t) and <fr for the action via the Dyson equa- 
tion. This procedure is repeated until convergence is 
reached. The selfconsistency equation for the condensate 
takes the simple form 



* = (b(r)) Simp . 



(32) 



The new hybridization function A(r) is obtained in 
the following way: From Eq. ([3~Tj) and Eq. (JS) we obtain 
the connected Green's function via 



G c (r) = G(t) + 



(33) 



We then Fourier transform G c (r) to obtain G c (iu> n ). In 
Appendix [C] we show how to obtain accurate Fourier 
transforms in spite of a finite number of measurement 
time-steps. With this we calculate the matrix selfenergy 
via the Dyson equation 



T,(iu n ) = G l (ioJ n ) - G c 1 (iuj n 



(34) 



Here G 1 (iw„) is the bare Green's function which is re- 
lated to the hybridization function A(iw n ) via 

A(ioj n ) = -iuj„a3 - jll + Gq (iu) n ). (35) 

The parameter fi = fi— (e) is chosen such that A(iw n ) — > 
in the limit u n — > oo. As we only consider symmetric 
density of states ((e) = 0) here, we write /x from now on. 
Eq. ([33)) allows us to rewrite the Dyson equation as 

S(iw„) = iu> n a 3 + ixl + A(ico n ) - Gj 1 (iw„). (36) 

Employing the DMFT approximation that the lattice 
selfenergy coincides with the impurity selfenergy, i.e. the 
selfenergy loses its momentum dependence: E(k, iw n ) — 



H(iu>n), we can calculate the fc-summed (or local) lattice 
Green's function with 



Gi a tt(^„) = ^ \iUn<T3 + (M ~ e k)l - H(iu r , 



. (37) 



where £k is the dispersion of the lattice. For some dis- 
persions £k it may be more convenient to transform the 
summation over wave vectors k into a integration over 
the density of states D(e). With D(e) = J2k^( e ~ £ k) 
Eq. (f3"T)l transforms to 

Gi a tt(iwn) = f deD(e)[iu) n <T3 + (fj, - e)I - 'E(iu) n )]~ 1 . 

(38) 

We now use the Dyson equation again to calculate the 
new updated hybridization function 

A(ioj n ) = -ioj n crs - A*I + S(iw n ) + G{J. t (iuj n ) (39) 

and the new value for k: 

K = zt- A u (iw„ = 0) - A 12 (iuj n = 0). (40) 

After an inverse Fourier transform we obtain A(r) and 
the selfconsistency loop is closed. We then solve the im- 
purity problem again with the updated action, until con- 
vergence is reached. 



A. Convergence 



In this section we 
cur in the iteration 
The first problem is 
the solution is still 
Hilbcrt transformat 
frequencies. The ori 
by writing Eq. (|38f 

Gi at t(«Wn) 



discuss two problems which may oc- 
process but can easily be overcome, 
that in the first few iterations, when 
far from the converged result, the 
ion may diverge for some Matsubara 
gin of the problem can be understood 
for the individual components 



deD(e) 



Gi a tt(iw„) = £ / deD(e) 



|C-e| 2 -£ 2 
1 



IC-e| : 



(41) 



where £ = iuj n + /i — E(«w„), E (S) are the diagonal (off- 
diagonal) components of the matrix selfenergy and we 
have used that E = E*. Obviously, the denominator of 
Eq. (|4"Tj) can vanish for certain values of E and E. If the 
divergence is due to the statistical errors on E and E, the 
problem can be cured by running more accurate simula- 
tions. However, another reason for a divergent integral 
can be that the initial approximations for A(r) and $ 
are unphysical ( "too far away" from the converged solu- 
tion). In almost all cases this problem can be avoided 
by choosing suitable initial hybridization functions A(r) 
and condensates 4>. We found that a good choice for the 
initial $ was the static mean-field value, corresponding to 
A(r) = 0. For this value of $ we then calculated G c (t) 
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and obtained our initial hybridization function A(r) us- 
ing the lattice Hilbert transform. Should the integral still 
be divergent, which only happens in the first few itera- 
tions, one can shift (increase) the value of £ such that 
the integral becomes well defined. Another strategy to 
fix this problem is to apply a damping procedure to the 
selfenergy in the first few iterations. 

The second problem is related to the convergence prop- 
erties of the condensate 3>. Sufficiently far away from 
the SF-Mott transition we typically need about 10 to 20 
iterations to reach convergence. Close to the SF-Mott 
transition, however, the convergence of 3? is very slow, 
and sometimes up to 500 iterations are needed to reach 
a converged solution. With a runtime of a few minutes 
per iteration this translates into a large numerical effort. 



Such slowing down of the convergence close to the phase 
transition can be overcome by using overrelaxation of 
the condensate The updating for $ is changed from 
Eq. © to 

* = a((b(r)) Simp - * old ) + * old , (42) 

where 3?oid is the condensate from the previous iteration 
and a > 1. 



V. MONTE CARLO METHOD 

To formulate our Monte Carlo Method we start by 
writing the impurity action Eq. (|29[) in non-matrix form: 



S imp = - / drdr' [6(r)F(r - r'^V) + b\r)K{r - t')&V) + b(r)K*{ T - r')b( T ')] 



dTn(r) + I 



(43) 



With this action we expand the partition function Z = Tr[Te s ™ p ] in powers of the hybridization functions F, K, 
K* and the source fields (j> and </>* . This leads to the series 

,/3 r fi ,/3 r P r fi r fi ,/3 ,/3 

Z = zZzZ ***■•••/ d <"A dr«... dr'< K dr>« . . ./ 

" ~ JO Jr£ , JO JO Jo Jt£ , Jo Jo 



n m 

rP 



x / dr{ 



K * 



K ' 



dr[ 



IK* 



IK" 



drf 



dr* / dr( 



dr* 



Teft J* Mr« ^KrJ-H^t^) . . . b{r^ F )b\r'i F )b{rr)Krr) ■ • ■ ) 

Xb\rf)b^rn . . . b\r^ K p{r^ K )b\rt) . . . b\r^)b{rf ) ■ ■ ■ &(<.)|n) 
xF(r[ - r[ F ) . . . F{t? f - r£)tf(r* - r[ K ) . . . K(t« k - r^J 

xK*(rf -r[ K *)...K*{T%* K „ - r£.) K m * +^^(rf ) . . . «^(<)0(rf * ) . . . 0(r^ ,), (44) 

where |n) denotes the eigenstate with n bosons and m = (mj?, fflx , m^* , )■ We can now sample individual 

diagrams with weight 



v ,( n . T F T F IF IF K K IK IK . K* K* IK* ,K' . d> A d>' T 4>' \ 

77 



Te ^ -mM-cr/a//" ^JWrj-il^) . . . b { r ^)b\r'Z)b{rf )■■■ K^m*. ) 

x ftt (7 f ) fc t (f f ) . . . b\r^ K )b\r^ K )b\rf) . . . b\r^)b(rf). . . &(<.)|n) 
x F(rf - Tf ) . • . F(r^ F - t£)* (rf - r[ K ) . . . K(t« k - r>« K ) 

x K* (rf -r[ K ')... K* (r,f K , - t'£, ) K m *^*' <f>(rf) . . . ^ )<j>{rf )... ^ ). 



(45) 



These contributions, illustrated in Fig. [TJ can be rep- 
resented by a collection of nif + 2niK<- + m<t> = nip + 
IrriK + m^. creation and annihilation operators on the 
imaginary time interval [0, /?). Hybridization functions F 



connect mp pairs of creation and annihilation operators, 
off-diagonal hybridization function K (K*) connect vhk 
(rax*) pairs of creation (annihilation) operators, while 
m<f> {jn$* ) creation (annihilation) operators are linked to 



source fields cj> (<fi*)- The configuration is fully specified 
by additionally giving the occupation n of the impurity 
at times r = 0, which in combination with the collection 
of operators determines n(r). 



F F 



<j>* <J> <j>* <J> 




FIG. 1: (Color online) Diagram corresponding to perturba- 
tion orders tuf = 1, rriK — 1, ttik* = 1, m^, = 2, m^,* = 2 
and n(r = 0) = 2. The hybridization function F determines 
the amplitude for transitions of bosons from the impurity into 
the normal reservoir, while operators coupling to source fields 
<j> and <j>* represent transitions between the impurity and the 
BEC reservoir. The off-diagonal hybridization functions K 
and K* , present only in the BEC phase, represent the am- 
plitudes for creating or annihilating two bosons at different 
times. 



1. <0|O |p> — < |«O|p> 



F F F 
2. a) (Q\^) 0|p> « ► <o|f-Q 0|p> 



2. b) <o| 



2. c) <0| 



K 



K 



K 



K 



|p> — <o|*-^>S-<> |p> 



K* K* 

|p> — <o| ^b^b |p> 

4> $ 



3. < | J=i |p> — <0h— B) |p> 



A. Updates and detailed balance 



F F F F 

4. < |£3> i^lp) - — ► <o|i-cT*-6|$ 



An ergodic sampling of all possible diagrams requires 
the following updates: 

1. insertion/removal of a pair 6(r)F(r — t')&'(t'), 

2. exchange of the bath type: 

a) b(T)F(T-T'p(T r ) o K( /)*6(t)k^(t'), 

b) b(T)K*(T-T')b{T') O K^*6(r)K0*6(r'), 

c) ^(t)K(t - T>t( T ') o «;0&t( T )«;06t( r '). 

In order to improve the efficiency additional moves 
can/should be used. Updates which may improve the 
efficiency are shifts of operator positions, moves which 
reconnect the hybridization lines and moves which in- 
crease/decrease n by one. For example, in the normal 
or Mott phase, where <j> = K = 0, or close to the SF 
transition, where <p is small, it is useful to have an ad- 
ditional move which reconnects two i^-lincs. The reason 
for this is that update 1 is the most time consuming and 
no F-lines can be reconnected via source fields <f>. Deep 
in the Mott phase the insertion of F-lines is strongly 
suppressed. In order to have an efficient sampling, we 
now need an update which increases/decreases the oc- 
cupation in the whole imaginary time interval without 
inserting new operators. A graphical representation of 
the updates is shown in Fig. [3J 

In order to satisfy detailed balance, we decompose the 
transition probability to go from state i to state j as 



p(i^j)=P woh (i 



(46) 



FIG. 2: (Color online) Illustration of the different updates 
described in Sec. IV Al 



The only move which changes the number of operators 
in the imaginary time interval, and therefore the total 
perturbation order to, is the move to insert or remove 
an F-line. If we choose to insert an operator pair (in- 
line) at random times r and r' drawn from a uniform 
distribution [0,/3) and propose to remove this pair with 
probability l/(mp + 1) the proposal probability becomes 



p prob (m F ^TO F + l) 
p prob (TO F + I — > nip) 



2dTdr' 
2 

TO p + I 



(47) 
(48) 



where p prob (p a 
this move. 



is the probability to propose (accept) Eq. (|45[) by WTr{^] T \ \ ■ 



The factor 2 comes from the fact that we have to 
decide if the occupation number is changed in the in- 
terval between the two operators (with length |r — t'|) 
or in the interval which winds around r = ft (with 
length /3 — \t — t'|). If we choose the interval which 
winds around r = (3 this will change the occupation 
n at time r = 0. Denoting the factor (n\ . . . \n) in 



), the 

detailed balance condition for inserting/removing a pair 



r' F 
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b{r)F{T -T')tf(r') becor 
p acc (mp —> 7Ji p + 1) 



1 



p acc (rriF + 1 — > ttif) 

I I F 
W T r(n -,T{ ,.. 



mp + 1 



F(r-r') 



w Tr (n]T{,.. 



T F T iF 

• ' raw ' 1 ' 



r'F 



•0 



(49) 



where n' and n are different only if the new _F-linc spans 
over time r = (and r = j3). 

For moves which increase/decrease the occupation of 
the impurity by one we have p prop = 1, so the acceptance 
probabilities satisfy 



p acc (n 



1) 



p acc (n + 1 — > n) 
WTr(n + 1; rf, . 



J ' rriF ' ' 



WTr(n;T{ 



T F rr-lF 



•) 



(50) 



Moves which interchange hybridization lines are easi- 
est, since the operator trace term WTr does not change. 
The probability to propose such a move is given by 
the inverse of the number of possibilities to choose the 
line (or lines) to exchange. For moves which exchange 
b(r)F(T - t'P(t') O K0*6(r)K06t(r') the detailed bal- 
ance condition reads 



p acc ({m F ,m^,m^} -s> {m F - l,m^ + l,ny« + 1}) 



p acc ({rriF — 1, ™0 + 1, 



rricf,* + 1} 

K(j)*K<j) 



{m^ + l)(m r +l)F(r-T') 



(51) 



For moves which change an off-diagonal hybridization 
line into two condensate lines we have to consider that 
the X-lines do not have a defined direction like the F- 
lines. Therefore the probability to choose one X-line in 
a state with ttik X-lines is given by pP rob = 2/rriK- Sim- 
ilarly, the probability to choose two </>-lines in a state 
with 7Ti0 0-lines is given by p prob = 2/{rn^(rn^ — 1)). 
The detailed balance condition for moves which exchange 
(t)K(t— r')6t(r') o K<pb^ (T)K<ptf (r') is therefore given 
by 



p acc ({m K ,m } 



{m K ~ l,m^ + 2}) 



(52) 



p acc ({m K - 1, + 2} -> {m K , n^}) 
m^- K(f>K(f) 
(m + 2)(m^ + l) A'(r-r')' 

and similarly for the complex conjugate. 

B. Measurement of observables 



By taking functional derivatives of the partition func- 
tion with respect to either F(t—t'), K(t—t'), K*(t—t'), 
or both 4>*(t) and 4>{t') one can calculate the diago- 
nal and off-diagonal parts of the Green's function matrix 
G(t), and the condensate order parameter </>. 



For the condensate the functional derivative yields 

where (A)mc means that the quantity A should be aver- 
aged over all configurations obtained in the Monte Carlo 
sampling. Due to the time independence of <fi this reduces 
to 



<f>=(b(r))s lmp = 



(54) 



One can see from Eq. (J54J) that one needs a condensate 
in order to measure a condensate. Therefore one always 
chooses (f> ^ as initial value in the simulation. 

The diagonal and off-diagonal Green's function can be 
measured in two different ways. Either by taking the 
functional derivative with respect to F(t — r'), K(t — t') 
and K*(t — t') or by differentiating with respect to both 
4>*{t) and 0(t'). This yields 



G(r) 



G(r) 



-(Tb(r)b\0)) 



Si, 



r /F^ 



-(Tb(r)b(O)) 



) / MC 



-K ' 



r rK* 



) / MC' 



(55) 



in terms of the hybridization and 



G(r) 



G(r) 



(T6(r)6t( )) Simp 

(EE {T - T ' ~ Tl ) 

i=i j=i 
-(T6(r)6(0)) s 

"*» rn 



^K<\)*K(\) I MC 



g 'g S(r,rf-rf ) 
i=i &i=i 



(3n(j)*K(f)* I mc 



(56) 



in terms of the condensate with 



6(t,t) 



5{t - f) for f > 0, 
5{t-t-P) forf <0, 



and similarly for the adjoint. The end point G(/3_) can 
be accurately measured through the density, 



'MC 



1 f 

drnir, 



MC 



(57) 



and G(0+) = G(/3_) — 1. In the condensate phase the ex- 
pansion order of 4> and <fi* is much higher than the expan- 
sion order of the hybridization, see Sec. IV Dl Therefore, 
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the Green's functions are measured by using the conden- 
sate. Close to the SF transition, where <j> is very small, 
it is more efficient to measure the diagonal Green's func- 
tion according to Eq. (|55|) . In the normal or Mott phase, 
where </> = 0, the diagonal Green's function can only 
be calculated in this way. In practice the off-diagonal 
Green's function is always measured with 0, even close 
to the SF transition. This is because the expansion order 
of K and K* is always lower than m^. 

The kinetic and potential energy can also be evaluated 
easily and very accurately. The potential energy is given 
by 



E 



U 



pot 



(n(n — 1)) — //(n), 



(58) 



and can be computed directly in the Monte Carlo simu- 
lation. As shown in Sec. MI Dl the kinetic energy in the 
DMFT approximation is given by -Ekin = ^\ —1 - Equa- 



tion ([24]) in the presence of a constant, homogeneous con- 
densate then yields 



E 



kin 



2^TrJ] [G c (w„)i?(G c K)) - I] - zt<f> 2 . (59) 



we obtain 



By using Eq. ([22]) and A 



Skin = ^Tr^ [G C K)A K)] - ztcf. (60) 

n 

Going back to imaginary time, we can write this in the 
individual components as 



(i.e., all diagrams obtained by linking these operators in 
different ways by hybridization functions) results in a de- 
terminant, and was the essential step which allowed to 
suppress the sign problem and to formulate an efficient 
algorithm. In a bosonic model, a similar summation of 
diagrams would lead to a matrix permanent whose eval- 
uation is ^P-complete. The only efficient way of evalu- 
ating permanents is their stochastic sampling, which is 
exactly our algorithm of individually sampling the differ- 
ent diagrams (as illustrated in Fig.Q} instead of summing 
them up expicitly. 

Since the off-diagonal hybridization function K is neg- 
ative, some diagrams have negative weights which leads 
to a sign problem. However, as we will show later, this 
becomes an issue only at very low temperatures in the 
presence of a condensate and does not prevent an accu- 
rate computation of phase diagrams and dynamical quan- 
tities. 



C. Solver test 

To check if the diagrammatic sampling and measure- 
ment procedure have been implemented correctly we 
compare the QMC-result with exact diagonalization. In 
a Hamiltonian formulation, the impurity model can be 
written as 

#im P = J2 [Viialb + aib^+Wii^b + ajb^ + J>a| 



-/in + ~^ n { n ~ 1) — K>{4>*b + <j>b^). 



(64) 



E kh 



<It[F(t )G c (t) + K(t)G* c (t) + K*(t)G c (t)} 



-ztcjr, 



(61) 



or in terms of the full Green's function 







Skin = / cIt[F(t)G(t) + K(t)G*(t) + K*(t)G(t)} 

>0 



-K(/) 2 



(62) 



By plugging Eq. (j54"| and Eq. (|55|) into the above expres- 
sion we see that the kinetic energy is directly related to 
the expansion order via 



E, 



(63) 



where m to t — mp + vtlk + tuk* + (w^ + m c f ) *)/2 is the 
total number of operator pairs in the interval [0,/3). 

This algorithm is a bosonic version of the hybridiza- 
tion expansion method of Ref. [l^. Aside from the con- 
densate term and the anomalous hybridization, the main 
difference lies in the fact that for bosons we don't use 
any summation of diagrams. In the fcrmionic case, the 
analytical summation of all diagrams corresponding to 
a given sequence of creation and annihilation operators 



The hybridization functions F and K in Eq. (|29[) are re- 
lated to the hybridization parameters Vi and Wi through 



F{iu n ) = J2 
i 



ViVi WiWi 



ei + iuj n ei - UJ r , 



ViWi vm 



(65) 
(66) 



For the test, we represent the bosonic bath by a finite 
number nbath of levels with creation (destruction) oper- 
ators a, (a;) and energies e;. By Fourier transforming 
Eq. (|66|) we get 



V, 2 



f M = £ - i 
k{t) I'gV W, 



ir, 



1 -e-^P 



2 ^ \ e^P - 1 1 - e~ e iP 

1=1 



(67) 



For the comparison we choose one orbital (nbath = 1) 
and the following parameters: t = 1, /3 = 1, /a = 20, 
U = 20, k = 6, <f> = 1, e = 1, V = 1 and W = -0.2. 
In Fig. [3] we show the diagonal and off-diagonal Green's 
function calculated with exact diagonalization (ED) and 
Monte Carlo simulations. 
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o Monte Carlo 

— Exact Diagonalization 





0.02 



0.025 



0.03 



0.035 
t/U 



0.04 



0.045 



0.05 



FIG. 3: (Color online) Comparison of the diagonal and off- 
diagonal Green's function calculated by Monte Carlo (black 
circles) and exact diagonalization (red curve) for a hybridiza- 
tion function which can be described by a bath with one or- 
bital. For parameters, see text. 



FIG. 4: (Color online) Mean perturbation order vs. t/U at 
filling n = 1 and inverse temperature fit = 2 for the 3d cu- 
bic lattice. The superfluid phase sets in where m</> and mis- 
become non-zero. 



The density and condensate order parameter are 
tied = 1-6762 and <^ed = 1-07976 for ED, and (umc) = 
1.67619(3) and (0 MC ) = 1.07974(9) for the Monte Carlo 
simulations. Note that for this set of parameters the av- 
erage sign (s) in the Monte Carlo simulation is (s) = 
0.40695(7). The perfect agreement with the ED result 
shows that the diagrammatic sampling and measurement 
procedure have been implemented correctly. Since the 
model we considered here contains both diagonal and 
off-diagonal hybridization terms and since a bath with 
a finite number of levels is as difficult to treat as any 
other bath, this serves as a nontrivial test for the Monte 
Carlo solver. 



D. Perturbation order and average sign 

In this section we show how the expansion order of the 
hybridization function scales with interaction and how 
the sign problem, caused by the off-diagonal hybridiza- 
tion K scales at the Mott insulator (or normal phase) to 
SF transition. Since we sample the operator configura- 
tion in the imaginary time interval [0, /3) the perturbation 
order grows roughly proportional to /3 in all phases. 

In Fig. 2] we show how the mean perturbation order 
grows as one goes from the normal phase to the SF phase 
at filling n = 1 . In the normal phase we have only contri- 
butions from the diagonal hybridization F and the per- 
turbation order tuf decreases with increasing U. In the 
SF phase the perturbation order grows rapidly with de- 
creasing interaction U /t, mainly because of the conden- 
sate contribution and the off-diagonal hybridization 
contribution vhk- 



A 
s= 
CO 




0.02 0.025 0.03 0.035 0.04 0.045 0.05 



t/U 

FIG. 5: (Color online) Average sign vs. t/U at filling n = 1 
for the 3d cubic lattice at different temperatures. The sign 
is negative (positive) for odd (even) expansion orders in vtik ■ 
Hence, where the average sign starts to deviate from 1 is in- 
dicative of the phase transition into the superfluid phase. 



Figure [5] illustrates the sign problem we encounter 
when going from the Mott phase to the SF phase. 
At the phase boundary where the condensate vanishes 
(and therefore K vanishes) the sign problem disappears. 
Therefore the sign problem does not prevent an accurate 
computation of phase diagrams. The sign problem is only 
severe deep in the condensate phase. 
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VI. SIMPLE LIMITS 
A. Non-interacting bosons 

For an ideal Bose gas (U = 0), the total number of 
particles is given by 



^ = I0| S 



rk- 



D(e) 



,0(e-M) _ I 



le < zt, 



(68) 



where zt is the half-bandwidth of the lattice. Just like 
for an ideal Bose gas in continuous space, the chemical 
potential fx has to be lower than the bottom of the band, 
\i < —zt in order to keep the number of particles finite. 

For temperatures T > T c , the condensate vanishes, 
1 <f>\ 2 = and the above equation determines the density 
n as a function of the chemical potential fx. For T < T c 
the chemical potential must be pinned at the lower edge 
of the band, /i = — zt in order to have condensation. The 
total number of particles is then a function of the con- 
densate density \(f>\ 2 . In the non- interacting limit the off- 
diagonal hybridization functions K and K* vanish and 
the B-DMFT equations equations become exact. The 
impurity action now takes the simple form 

rP rP 
S imp = - / dTdT'b{T)F{T -t')^{t')- / drn(r) 
Jo Jo 

(69) 







-K / dr[(j>* {r)b{T) + cj)^ {t)\. 
Jo 



This is a quadratic action that can be solved analyti- 
cally. The solution for the non-interacting Green's func- 
tion and hybridization function is given by 



G(iuj n ) = G (iu} n ) = / de- 



Die) 



li + G 



iu) n + [i — e 
-iw n ). 



(M < **) 
(70) 



B. Static mean-field 



One obtains a selfconsistent mean-field theory (the de- 
coupling approximation) by substituting 

b\h = (bt)bi + bt(bi)-(b\)(bj) 

= 0(&t + 6j .)_02 (71) 

into our Hamiltonian defined by Eq. (TT]). If one drops 
the term — </> 2 , which is just a constant shift in energy 
one obtains the following Hamiltonian 

U 



this just gives the half-bandwidth k — 6t. This Hamilto- 
nian can be expressed as a matrix in the occupation num- 
ber basis (truncated at some maximum occupation) and 
solved by exact diagonalization. One chooses an initial 
value for the condensate <f> and determines 4> iteratively 
by solving cf> = (b) until convergence is reached. 

By using B-DMFT we can reproduce the static mean- 
field results by setting the hybridization function to zero, 
i.e. A(t) = 0. Since there is no hybridization Eq. (|30|) 
reduces to k = zt. 



VII. NUMERICAL RESULTS 



In this section we present results for the Bose-Hubbard 
model on a 3d cubic and 2d square lattice obtained with 
B-DMFT. All our results are compared to other methods 
like static mean-field theory, worm-type quantum Monte 
Carlo (QMC) simulation on a lattice of up to 40 3 sites^^ 
and with a recently developed numerically exact method 
on the Bethe latticed Since the Bethc lattice can be also 
directly simulated with B-DMFT we will show a direct 
comparison between B-DMFT and the numerically exact 
solution. 




oooooooooooooo 
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FIG. 6: (Color online) Kinetic energy of the Bose-Hubbard 
model on the 3d cubic lattice as a function of temperature 
for fi/t — 8 and U/t = 20 (n w 1). Results obtained from 
B-DMFT (black circles) are compared to lattice QMC (blue 
diamonds) and to static mean-field theory (black dashed line) . 
Inset: Energy difference from the QMC data for the same 
parameters. The QMC results are obtained on a lattice with 



Tj ji Wit u \ i i -\\ 10 sites except close to the transition (4 < Tit < 6) where 

— ' 1 9 ^ — ' ^ — ' 40 sites were used. Krror bars are smaller than the svmbol 

size. 



(72) 

where n = -\ t = zt is the hopping term summed over 
the nearest neighbors. For the 3d cubic lattice (z = 6) 



Error bars are smaller than the symbol 
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FIG. 7: (Color online) Total energy of the Bose-Hubbard 
model on the 3d cubic lattice as a function of temperature 
for fi/t = 8 and U/t = 20 {n ~ 1). Results obtained from 
B-DMFT (black circles) are compared to lattice QMC (blue 
diamonds) and to static mean-field theory (black dashed line). 
Inset: Energy difference from the QMC data for the same pa- 
rameters. The QMC results are obtained on a lattice with 
10 3 sites except close to the transition (4 < T/t < 6) where 
40 3 sites were used. Error bars are smaller than the symbol 
size. 



A. 3d cubic lattice 

1. Kinetic and total energy 

In this section we present the results for the 3d cu- 
bic lattice. Since DMFT can be considered as an ap- 
proximation to the kinetic energy functional, Eq. (|25p . 
it is interesting to see how the kinetic energy per site 
obtained from B-DMFT compares to results from lattice 
Monte Carlo simulations^ and static mean-field theory. 
In Fig. [S] we show the kinetic energy as a function of tem- 
perature as one goes from the SF to the normal phase. 
In the case of static mean-field theory, the kinetic en- 
ergy is just given by Skin = —zt(f) 2 , where cj) = (b) is 
the condensate order parameter. In the ground state 
regime this gives a good approximation of the kinetic 
energy. For T > T c , where = 0, the kinetic energy van- 
ishes since hopping of the normal bosons is completely 
neglected. The agreement of the B-DMFT results with 
the exact QMC results is excellent over the whole tem- 
perature range. Only close to T c there is some small 
deviation. For the QMC simulation we used 10 3 lattice 
sites, except for temperatures in the range 4 < T/t < 6 
where 40 3 sites where used. In Fig. [7] we show the to- 
tal energy as a function of temperature. The remarkable 
accuracy of the B-DMFT result for the total energy at 
all temperatures implies that entropies may also be com- 
puted reliably using B-DMFT. 

In Fig. [5] we show the finite temperature phase dia- 



gram (top panel) and the ground state phase diagram 
(bottom panel) for the first and second Mott lobe of the 
Bose-Hubbard model on a 3d cubic lattice and compare 
results obtained with B-DMFT to exact results from lat- 
tice QMC simulations, the exact solution for the Bethe 
lattice with coordination number z = 6^ and to static 
mean-field results^ The SF phase is characterized by a 
finite value of <f> = (b), while we have <f> = (b) = in the 
Mott insulating and normal phase. For the calculation of 
the ground state phase diagram we used fit = 2, which 
is shown in Fig. [5] to be a sufficiently low temperature. 
Based on simulations done at (3t = 4 and (3t = 8 we found 
that any systematic error is smaller than the statistical 
error. The excellent agreement between our B-DMFT 
results and the full solution of the Bose-Hubbard model 
shows that the Mott-transition is a local phenomenon, 
well described by a momentum-independent sclfcnergy 
and that the condensed bosons are accurately described 
by a uniform condensate. 



2. Critical exponents 

We now examine the critical exponent of the order 
parameter (f>. In Ref. [10 fermionic DMFT was em- 
ployed to study the superfluid state. The temperature 
dependence of the condensate order parameter A goes as 
A ~ |(T C - T)/T c f in the vicinity of the critical tem- 
perature T c , with ,3 = 1/2 being the mean-field expo- 
nent. The symmetry breaking happens on the mean- 
field level in the two-particle channel given by a spin-up 
and a spin-down fermion. For any fermionic operator 
we have that (c) = such that condensation can occur 
the earliest in the two-particle channel, e.g., (c|c^) or 
(cj-cj,) may be non-zero. For a bosonic operator, a sin- 
gle operator b may already have a non-zero expectation 
value, (b) ^ 0, leading to condensation at the one-body 
level. The B-DMFT action allows for this, while at the 
two-body level (i.e., the selfconsistcntly determined hy- 
bridization terms) the one-loop correction to the conden- 
sate will change the critical exponents from their mean- 
field values. The exponents we obtain are therefore not 
universal but depend on the parameters i, /i, U , (3 and 
the lattice dispersion ek- Whenever k — const, we recover 
the mean-field exponents. This happens of course in the 
static mean- field limit (k = zt) and for non- interacting 
bosons (« = — Gq 1 (iu) n = 0)), where the mean-field ex- 
ponent p = 1/2 is exact. In Ref. [49l mean-field exponents 
were also found for the exact solution on the Bcthc lat- 
tice. 

In Fig. [9] (main panel) we show the condensate order 
parameter as a function of temperature obtained by 
B-DMFT and compare the results to lattice QMC and 
static mean-field theory. Close to the critical point the 
condensate goes as ~ |(T C - T)/T c |' 3 for T < T c . The 
QMC results are obtained on a lattice with 40 3 sites, 
which still leads to a rounding of the data compared to 
the thermodynamic limit. They are based on the k = 
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FIG. 9: (Color online) Main panel: Condensate order pa- 
rameter 4> as a function of temperature of the Bose-Hubbard 
model on the 3d cubic lattice for fj,/t — 8 and U/t = 20 
(n w 1). Results obtained from B-DMFT (black circles) 
are compared to lattice QMC (blue diamonds) and to static 
mean-field theory (black dashed line). Inset: Zooming in on 
the region close to the SF-normal transition showing the crit- 
ical exponent ft. Here 4> is plotted as a function of T c — T 
for T < T c . The full lines are fits to <j> = A\T C - Tf (B- 
DMFT: p = 0.194(2), T c /t = 4.365(3), At = 0.6463, QMC: 
P = 0.35, T c /t = 4.43(3), At = 0.6517, Mean-field: p = 0.5, 
T c /t = 6.5661, At = 0.535). The QMC results are obtained 
on a lattice with 40 3 . Error bars are smaller than the symbol 
size. 



FIG. 8: (Color online) Top panel: Phase diagram (superfluid 
to normal liquid transition) of the cubic lattice Bose-Hubbard 
model in the space of interaction and temperature for n = 1. 
The dashed line shows the static mean-field result, the red 
curve the exact solution for a Bethe lattice with coordination 
number z — 6 (Ref. |49T ) and the blue curve with open dia- 
monds the QMC result from lattice simulations (Ref. l48l) . The 
black line with open circles corresponds to the B-DMFT so- 
lution, which yields a second order transition. Bottom panel: 
ground-state phase diagram in the space of t/U and (J./U, 
showing the first two Mott lobes surrounded by superfluid. 
The B-DMFT phase boundary was computed at fit — 2. Er- 
ror bars are much smaller than the symbol size. 



3. The weakly interacting Bose gas regime 

When bosonic field theories are expanded in U (cf. 
the weakly interacting Bose gas of Ref. [44]), the effect 
of the chemical potential is non-perturbative. The chem- 
ical potential is negative for the ideal Bose gas, and has 
to change sign in the presence of repulsive interactions. 
There is an implicit relation between the condensate den- 
sity no = \(j)\ 2 and the chemical potential, which follows 
from the condition of thermodynamic equilibrium,— 

\dn(T = 0,V,»,n )] _ n 



component of the Fourier transform of the equal-time 
density matrix. By plotting ^ as a function of T c — T 
in a log-log plot we can extract the critical exponent (3 
from the slope of the line as T — > T c , which is shown 
in the inset of Fig. O Static mean-field theory of course 
gives the mean-field exponent (3 = 1/2. The exact model 
belongs to the 3d XY universality class with the exponent 
13 w 0.35ri while the exponent obtained from B-DMFT 
for these parameters is f3 ss 0.19. 



leading to /i = (0| |0), with H mt the interacting two- 
body terms and |0) the ground state (finite temperature 
extensions also exist^). In any expansion order this re- 
mains valid and can be worked out to yield 

fx = S(fc = 0,w = 0) - E(fc = 0,u> = 0), (74) 

which is the famous relation for a gaplcss spectrum first 
derived by Hugenholtz and Pines^ 

The Hugcnholtz-Pincs relation does not hold in our B- 
DMFT formalism. To start the discussion, let us consider 
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FIG. 10: Relative deviation from the Hugenholtz-Pines rela- 
tion, (5/E(w n = 0) with 8 = E(cj n = 0) — E(cj„ = 0)— /x— zt, of 
the Bose-Hubbard model on the 3d cubic lattice as a function 
of temperature. The figure shows the transition from the SF 
to the normal phase for U/t — 10 at filling n m 1. Error bars 
are of the order of the symbol size. 



a mean-field approach in which the condensate is deter- 
mined sclfconsistcntly but the hybridization functions are 
not. Then it is easy to see that the Hugenholtz-Pines re- 
lation is shifted to /i = S(cj = 0) — £(<J = 0) — zt. If, 
however, the hybridization functions are also determined 
sclfconsistcntly (as in the full B-DMFT scheme), then 
they will also depend on and influence the condensate 
density (and contain fi in them), leading to the break- 
down of the Hugenholtz-Pines relation. Indeed, although 
Ref . uses a skeleton approach for the description of the 
WIBG, the discussion of the Hugenholtz-Pines relation is 
provided in terms of bare Green functions in order not 
to mix the expansion orders. 

However, it turns out that the deviations from the 
(shifted) Hugenholtz-Pines theorem are minimal. This 
is shown in Fig. [TU] where we plot the relative deviation 
from the Hugenholtz-Pines relation at a value of U/t = 10 
across the SF to normal transition. One clearly sees that 
the deviation grows as one moves deeper into the SF 
phase. In the Mott and normal phases, where = 0, 
the Hugenholtz-Pines relation is not valid. 

As already mentioned in the introduction, it was shown 
in Refs.EHi that gaplcss field theories have selfcnergies 
that should be momentum dependent and, in particular, 
that the anomalous selfenergy has to go to zero for small 
momenta at zero frequency. In Beliaev's diagrammatic 
approach of the weakly interacting Bose gas, the first or- 
der selfenergies are momentum and frequency indepen- 
dent, and this condition is not satisfied, questioning Be- 
liaev's approach. Similarly, in the B-DMFT scheme the 
selfcnergies are momentum independent. In Fig. 1111 one 
sees that the anomalous selfenergy is non-zero for zero 
frequency, as could be expected from the explicit sym- 



FIG. 11: (Color online) Selfenergy of the Bose-Hubbard 
model on the 3d cubic lattice as a function of Matsubara 
frequencies in the SF phase, j3t = 1, fj,/t = 8 and U/t = 20 
((n) ~ 1). The curve with blue circles (diamonds) shows 
the real (imaginary) part of the selfenergy E(iw n ) and the 
curve with red triangles shows the off-diagonal selfenergy 
E(jWn). The dashed lines show the analytically known high 
frequency tail where E(iw„) = Eo + Ei/«o; n + 0(l/(icj„) 2 ) 
and E(«w„) = En + 0{l/{iu n ) 2 ) with E = 2U(n), Ei = 
U 2 (3{n 2 ) - (n) - 4(n) 2 + (6b) 2 ) and E = U{bb). 



metry breaking terms in the B-DMFT action. On the 
other hand, the arguments of Ref. HH for the weakly in- 
teracting Bose gas also hold for B-DMFT. In particular, 
the authors of Ref. argued that the system reaches 
its thermodynamic limit for properties such as the en- 
ergy and the entropy for small system sizes. In Figs. [51 E] 
we have already seen that the (kinetic) energy is repro- 
duced remarkably well over the entire temperature range 
in 3d, except in the fluctuation region near the normal- 
supcrfluid phase transition point. 

From the viewpoint that DMFT and B-DMFT sum up 
all skeleton diagrams built with local propagators only 
(which is an (asymptotic) expansion in U), we expect 
that our B-DMFT theory should recover the physics of 
the weakly interacting Bose gas (Bogoliubov Hamilto- 
nian), which is known to be an excellent approximation 
at low values of nU . It is the only limit for which it 
is not obvious that the B-DMFT formalism will be suc- 
cessful, while the other regimes of strong U or high tem- 
perature are well captured by a single-site action. In 
the skeleton approach to the WIBG theory^ the self- 
energies are to leading and subleading order given by 
E = 2(n)U and £ = (bb)U (see Fig. [TSJ. To this 
order of accuracy, the chemical potential is given by 
(jp) = 2{n)U— {bb)U rj {tiq)U. The latter approximation 
is made to keep the equations simple (it is not fundamen- 
tal) since we don't want to go into technical details of the 
WIBG theory here^ but it allows us to simply relate the 
condensate (no) to the chemical potential jjS 1 ^ ~ /i in 
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FIG. 12: (Color online) Selfenergy of the Bose-Hubbard 
model on the 3d cubic lattice as a function of Matsubara fre- 
quencies in the SF phase, fit = 4, fi/t = —5.75 and U/t = 0.5 
~ 0.5). The parameters are in the regime where the 
weakly interacting Bose gas (WIBG) theory applies, which 
predicts E (1) = 2{n)U and E (1) = (bb)U. The selfenergies are 
almost frequency independent and are in very good agreement 
with the WIBG predictions. 



FIG. 13: (Color online) Shown is the chemical potential as a 
function of U/t such that the density is approximately (n) — 
0.5 for the case of the weakly interacting Bose gas (WIBG), 
B-DMFT, quantum Monte Carlo worm-type simulations on 
a lattice of size L = 10 3 , and static mean-field theory. The 
temperature is T/t = 0.25. For low U, both static mean-field 
theory and B-DMFT approach the weakly interacting Bose 
gas theory well, while in the intermediate regime B-DMFT 
is doing a superior job for the chemical potential despite the 
local approximation for the self energy. 



this limit, after also taking the shift — 6i from the band 
edge into account. Because the chemical potential in B- 
DMFT is an input parameter and fully renormalized from 
the start, this is another way of understanding the devi- 
ations from the Hugcnholtz-Pines relation. 

We now discuss the numerical results for the weakly in- 
teracting Bose gas limit. In order to see some depiction 
effects, we choose parameters such that (n) = 0.5 and 
j3t = 4. The main result can be seen in Fig. [T^J where 
one sees that the selfenergy and the anomalous selfen- 
ergy become frequency independent and approach their 
WIBG values, which are precisely the leading terms in 
the high-frequency expansions used in B-DMFT: For the 
normal selfenergy it is S = 2(n)U, while for the anoma- 
lous selfenergy it is U (bb) . 

In Fig.[13]thc chemical potential is calculated such that 
the density is (n) = 0.5. The chemical potentials found 
by B-DMFT are in excellent agreement with the ones 
found by quantum Monte Carlo worm-type simulations. 
This shows again that the deviation from the Hugenholtz- 
Pines relation in B-DMFT is minimal, and goes to zero 
as U/t —> 0. We also compare the condensate density 
for these theories in Fig. 1141 a more challenging quan- 
tity than the total density. The agreement with quan- 
tum Monte Carlo simulations and the agreement with 
the WIBG theory for U/t < 0.5 is remarkable, and the 
difference with static mean-field theory is notable. Tech- 
nically, calculations in this limit are hampered by the low 
values of the selfenergies and the proximity of the band 
edge due to which very precise calculations are needed. 
For the density (n) =0.5 the sign problem is not the lim- 



iting factor (the average sign is approximately (s) = 0.7), 
but the situation deteriorates quickly for higher densities 
at these temperatures. 

To conclude this section, we briefly comment on the 
difference between the anomalous selfenergy found here 
(E^ = (bb)U) and the one found in Ref. (£W = 
(no)U) for a system in continuous space. The difference 
between the two is the diagram given by the convolu- 
tion of the anomalous Green function with the interac- 
tion line, at finite momenta (that is the second diagram 
for S^ 1 ) in Fig. [T5|l . In Ref. this diagram was not 
present in first order for the anomalous selfenergy since 
it was argued that interactions can to leading order be 
replaced by their k = values (Eq. (3.44), for which 
this diagram is zero. The diagram docs contribute to the 
chemical potential however for non-zero momenta, see 
Fig. 9 on page 22 of Ref. 0. In case of B-DMFT how- 
ever, we cannot make the distinction that the strength 
of the interaction is different for small or large momenta. 
Hence, to leading order there is still a contribution of the 
anomalous Green function to the anomalous self energy; 
in other words U(bb) is the correct value for the selfenergy 
in the WIBG limit for B-DMFT to leading order. 



B. 2d square lattice 

In this section we present results for the 2d square lat- 
tice. As in the previous section data from B-DMFT are 
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FIG. 14: (Color online) Shown is the condensate fraction as a 
function oiU/t such that the density is approximately (n) = 
0.5 for the case of the weakly interacting Bose gas (WIBG), 
B-DMFT and static mean-field theory. The temperature is 
T/t = 0.25. For low and intermediate U, B-DMFT is in 
better agreement with the weakly interacting Bose gas than 
static mean-field theory. 
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FIG. 15: First order skeleton diagrams of the selfenergy valid 
in the weakly interacting Bose gas regime, with results E 1 - 1 ' = 
2(n )U + 2G c (T = 0~)U = 2(n)U and E (1) = {n )U + G c {t = 
0)U = (bb)U. Note that in single-site B-DMFT all momenta 
are equal, and our interactions are instantaneous. 
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FIG. 16: (Color online) Top panel: Phase diagram of the 
square lattice Bose-Hubbard model in the space of interaction 
and temperature for n = 1. Bottom panel: ground-state 
phase diagram in the space of t/U and fi/U, showing the 
first two Mott lobes surrounded by superfiuid. The B-DMFT 
phase boundary was computed at fit = 8. The dashed line 
shows the static mean-field result, the red curve the exact 
solution for a Bethe lattice with coordination number z = 4 
(Ref. l49h and the blue curve with open diamonds the QMC 
result from lattice simulations fRef. [54ri. The black line with 
open circles corresponds to the B-DMFT solution. Error bars 
are much smaller than the symbol size. 



compared to lattice QMCj« static mean-field theory and 
an exact solution of this model on a Bethe lattice with 
connectivity z = 4.— In Fig. [IBlwe show the finite tem- 
perature phase diagram (top panel) and the ground state 
phase diagram (bottom panel) for the first and second 
lobe on the 2d square lattice. The ground state phase 
diagram was calculated at (3t = 8 which is sufficiently 
low as can be seen from the top panel of Fig. [THJ As 
one might expect, the agreement of B-DMFT with exact 
Monte Carlo lattice simulations is not as good in 2d as 
in 3d. We note here that in Fig. [16] we compare the con- 
densate transition obtained by B-DMFT to the SF tran- 



sition obtained by Monte Carlo lattice simulations. Since 
B-DMFT is just a dynamical extension to static mean- 
field theory we obtain (as in the static case) a finite con- 
densate (ft which is not present in a 2d system except at 
T = according to the Hohenberg-Mermin- Wagner theo- 
rem. Similar as in the discussion of Ref. |44J on the weakly- 
interacting Bose gas theory, Beliaev's diagrammatic tech- 
nique is useful in computing thermodynamic observables 
since they converge on short-range distances. The long- 
range fluctuations of the phase that ultimately destroy 
the condensate can be described at long wavelengths with 
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a separate hydrodynamical action. The diagrammatic 
technique simplifies hence enormously by working with 
explicit symmetry breaking and a genuine condensate in 
the calculation of thermodynamic observables4^ These 
ideas are illustrated for the kinetic energy in Fig. [T7l 
where at high and low temperatures the agreement with 
quantum Monte Carlo worm-type simulations is excel- 
lent, in contrast to static-mean field theory. Only in the 
vicinity of the transition point relatively small deviations 
are found. As in 3d, B-DMFT predicts a kink in the ki- 
netic energy curve at the transition point, and is hence 
not well equipped to describe Kosterlitz-Thouless phase 
transitions. Still, it gives a substantially improved finite 
temperature phase diagram compared to static mean- 
field theory (top panel of Fig. [16]) and a ground state 
phase diagram (bottom panel) in good agreement with 
the QMC result, 
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FIG. 18: (Color online) Density n and condensate order pa- 
rameter cf> = (b) as a function of fi/U for t/U = 0.02 and 
fit = 1 on the Bethe lattice with connectivity z = 4. The 
open black circles are the result obtained with B-DMFT and 
the solid red line the numerically exact solution (Ref. |4!|) . 
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FIG. 17: (Color online) Kinetic energy of the Bose-Hubbard 
model on the 2d square lattice as a function of temperature 
for n/t — 3.52 and U/t = 10 (n ft* 1). Results obtained from 
B-DMFT (black circles) are compared to lattice QMC (blue 
diamonds) and to static mean-field theory (black dashed line). 
Inset: Energy difference from the the QMC data for the same 
parameters. The QMC results are obtained on a lattice with 
100 2 sites. Error bars are smaller than the symbol size. 



C. Bethe lattice 



1. Finite connectivity 



tivity z the density of states is given by 

D(e) = ^S(e + zt)+D c (e), (75) 

with V the volume of the system. The continuum part 
D r reads 



Dr. 



^/4(z-l)t 2 -e 2 
2ir(zt 2 - e 2 /z) 



|e| < 2ty r z~l, (76) 



where t is the hopping between the sites. The isolated 
state at e = — zt matters for condensation and prevents 
the existence of Goldstonc modes in the symmetry broken 
phase.— In Fig. [18] we show a comparison of the density n 
and the condensate order parameter <f> between B-DMFT 
and the exact solution for a Bethe lattice with connec- 
tivity z = 4. As one can see in Fig. [T9] the agreement 
is excellent except for the condensate <j> very close to the 
SF to Mott transition. This comes from the fact that 
B-DMFT does not reproduce the mean-field exponents 
like the exact solution^ but has a smaller exponent in 
this case, i.e. the transition is much sharper. This ex- 
cellent agreement of the B-DMFT condensate <j> with the 
exact solution also shows that the definition of the con- 
densate used both in this work and in our previous Letter, 
Ref. [H, is correct, unlike the condensate defined in the 
Comment^ on our previous workp& which corresponds 
to the blue diamonds in Fig. Q1J] 



In this section we compare the solution of the B-DMFT 
equation on the Bethe lattice with the exact solution by 
Semerjian and co-workers4^ Since both methods are for- 
mulated in the thermodynamic limit this allows us to 
directly compare results and to check the accuracy of the 
B-DMFT formalism. For the Bethe lattice with conncc- 



2. Limit of infinite connectivity 

For systems without symmetry breaking, the isolated 
state in the non-interacting density of states is negligi- 
ble. It is customary in (fcrmionic) DMFT to rescale the 
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FIG. 19: (Color online) Condensate order parameter <f> = (b) 
as a function of fi/U for t/U = 0.02 and f3t = 1 on the Bethe 
lattice with connectivity z = 4 (blow-up of <f> in Fig. ll8l around 
fi/U = 2). The open black circles are the result obtained with 
B-DMFT, the solid red line the numerically exact solution 
(Ref. |49T) and the blue diamonds are results for the condensate 
<j> defined as in Ref. [5^. 



hopping in the limit z — > 00 as t = t/y/z in order for the 
density of states to remain finite. Equation (f76|) reduces 
then to the semi-circular density of states given by 



D(e) = - e\ \e\<2l (77) 

For this density of states the selfconsistency equations 
simplify since the Hilbert-transformation can be per- 
formed analytically. This allows us to formulate the 
whole selfconsistency loop in r-space, i.e. we do not have 
to go to Matsubara frequencies. We obtain a simple rela- 
tion between the hybridization and the Green's function 
given by 



A(r) = -PG c (t). (78) 

For bosons, however, the isolated state in the non- 
interacting density of states is not negligible. As soon 
as z > 2 the condensation temperature is non-zero, and 
one has then to keep track of the uniform mode^ Inte- 
ger scaling t = t/z is then the only possibility, and the 
static mean-field theory^ becomes exact for z — > 00 (for 
all temperatures). This integer scaling is also compatible 
with Eg. (f30|) such that the hybridization terms become 
identically zero for z — > 00, again in line with Refs. [i"ll55l 
but in contrast with the different scalings postulated in 
Ref. |M 



VIII. CONCLUSIONS 

In summary, we have derived the B-DMFT formal- 
ism in a number of different ways: (i) by approximat- 
ing the kinetic energy functional, (ii) by considering a 
cavity method or 1/z expansion, (iii) by using an effec- 
tive medium approach. We highlighted similarities and 
differences with previous approaches.— ~ 35 ' ' 49 We dis- 
cussed the Monte Carlo impurity solver in great detail 
including detailed balance of the updates, technical dif- 
ficulties one may encounter, and the measurement pro- 
cedure. We have shown results for the 3d Bose-Hubbard 
model, where the zero temperature phase diagram and 
the finite temperature phase diagram are in excellent 
agreement with quantum Monte Carlo worm algorithm 
results. We compared with the theory of the weakly in- 
teracting Bose gas in the regime of weak repulsion. The 
kinetic energy is reproduced at the 1 — 2% level over the 
full range of temperatures, except in the very close vicin- 
ity to the superfluid-to-normal transition where devia- 
tions are larger. We studied the critical exponents of this 
transition on the supcrfluid side, and found non-universal 
values depending on the parameters of the Bose-Hubbard 
model that are different from the mean-field values and 
from the universal ones found in the 3d XY universal- 
ity class. The reason is that the contribution of the hy- 
bridization terms can be understood as a one-loop cor- 
rection to the condensate. In two dimensions, the ground 
state phase diagram is still in good agreement with quan- 
tum Monte Carlo worm simulations, but at finite tem- 
perature the superfluid-to-normal transition is less ac- 
curately reproduced, because of the Kosterlitz-Thouless 
nature of the phase transition: in B-DMFT there is still 
a genuine condensate at finite temperature in 2d, and 
we miss the long-range fluctuations of the phase. How- 
ever, quantities that are local such as the energy are still 
remarkably accurately reproduced over the entire param- 
eter regime, just as in 3d. We also compared B-DMFT 
with the exact solution of Ref. on the Bethe lattice, 
and found again good agreement except in the very close 
vicinity of the phase transition points. The agreement 
over the whole parameter regime (except in the vicinity 
of phase transitions) was argued to be the result of fast 
diagrammatic convergence of thermodynamic quantities 
with system size, similarly as for the weakly interact- 
ing Bose gas system.— Our approach constitutes thus a 
natural extension of the Bogoliubov mean-field theory of 
the weakly interacting system (which is unable to find a 
Mott insulator) to a dynamical mean-field theory success- 
ful over the entire parameter regime of the Bose-Hubbard 
model. It shows that, when technical difficulties such as 
the asymptotic behavior of the series are properly taken 
care of or sidestepped, elusive diagrammatic expansions 
for bosonic systems are in fact promising, in line with 
the spirit of (bold) diagrammatic Monte Carlo^i studies, 
whether or not in combination with DMFT.— 

The virtue of developing a DMFT formalism for the 
Bose-Hubbard model lies in the possible extension of 
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the formalism to models where known numerical meth- 
ods fail. Prime examples are Bose-Fermi mixtures in 
two or three dimensions,— ~— where previous calcula- 
tions treated the spinless fermions within DMFT and the 
bosons at the mean-field levels The formalism can also 
be generalized to mixtures with spinful fermions. Clus- 
ter expansions are possible along the lines of Ref. |43|, but 
it remains to be seen whether the sign problem remains 
tolerable in the condensed phase. Developing the formal- 
ism for real-time applications can be considered along 
the same lines as for fermionic DMFT,— but here, too, 
it remains to be seen how the sign (as well as a complex 
condensate) behaves. 
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Appendix A: Quantum Cavity Method (1/z 
expansion) 



In a recent paper^ Semerjian et al. derived an ex- 
act solution for the Bose-Hubbard model on the Bcthc 
lattice (lattices with tree structure). Their final result 
is a selfconsistency equation that can easily be solved 
numerically. The paper mentions that B-DMFT can be 
considered as the first order approximation to this equa- 
tion. This is indeed the case, and leads to an equivalent 
derivation of the B-DMFT equations. The expressions in 
Ref. Hil are however for the cavity condensate and Green's 
function, not for the true (observable) condensate and 
Green's function. In this paragraph we write down their 
equations for the observable quantities and compare with 
our B-DMFT action^ 

Computations on a tree can benefit from the recursive 
structure of the tree. Following Ref. HH, the quantity 
Zi^j({n}) for two adjacent sites i and j is defined as the 
partial partition function for the subtree rooted at i and 
excluding the branch to j for given quantum numbers 
{n}. We note that a Suzuki- Trotter decomposition for 
the quantum action preserves the structure of the lattice 
for every time slice. Normalized quantities which can 
be interpreted as probabilities can then be defined as 
V^AM) = WW)/ £«„'}) ^WKD- We refer 
to Ref. 49 for more details. 

When taking an infinite number of Suzuki- Trotter 
slices, the cavity field on a Bcthe lattice of connectivity 
z can selfconsistently be expressed as^ 



Vc,v(b*,b) = ^-w(b*,b) /n»K.ii]'fc w WA)fflP t /V{6*(t)|>(t)+&(t)V&*(t) 
Zcav J i=i L Jo V i=i i=i 



with the on-site weight of the path given by 



)(b*,b) = cxp 



&*( T )(0 r - M )&(T) + ^n(r)(n(r)-l) 



(Al) 



(A2) 



r 



In Nambu notation this reads 
1 



7?cav(b) 



Z r 



-w(h)e 



(z-i)r(tb) 



= In 



rP 

2?b77 cav (b)exp / drb t (r)*(r) 
Jo 



(A3) 



This selfconsistency equation can only be solved exactly 
for U — 0. For the interacting case, U > 0, a solution 
can be looked for in the large connectivity limit, where 
in the leading order 1/z B-DMFT is recovered. This 



follows from an expansion of the generating functional of 
connected correlation functions T, 



r(#) 



rfr(b T 



#(t) 

^ P drdr'*t (T ) G ^ av (r-r')*(r') 



(A4) 
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The connected part of the two-point correlator satisfies 



bt(r)b(r'))cav-(b t (T)) cav (b(r')) c 



Gcav( T — T 

We now plug Eq. (|A4[) into Eq. (|A"3l) and obtain 
1 
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rP 
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G, 



(5 T a 3 - //I)5(t - r) 
-t 2 (z-l)G c c av (r-r' 



(A6) 



In the limit of large connectivity, one has to scale the 
hopping as t ~ 1/z, hence zt ~ 1 and £ 2 z ~ 1/z. 

We now proceed to the equations for the true fields, 
which cannot be found in Ref. |49U^ To this end, the true 
marginal rj(b) = ^w(b)e zr ^ tb ^ is needed. It is still ex- 
pressed in terms of the cavity field (and not the true field) 
and is obtained by replacing z — 1 by z in Eq. (|A6I) . We 
try to find now suitable expressions for the true connected 
Green's function and condensate, valid for the impurity 
problem 'imp' instead of for the cavity problem 'cav'. 
The action for the impurity problem is valid up to order 
1/z. Since the prefactor of the cavity Green's function 
is already of 0(1/ z), we can identify the impurity Green 
function with the cavity Green function without loss of 
accuracy, 

Gf mp (r-r') = G c cm (r-r') 

= (bt(r)b(r')) im p - <bt(r)) imp (b(r')> imp , 

(A7) 

where the average can now be taken with respect to the 
impurity action. 

The prefactor of the condensate in Eq. (|A6[) is of order 
unity, hence we need the impurity condensate to 0(l/z). 
We have 



r P 

Simp = Scav-* / dr*t b ( T ). 
Jo 



(A8) 



The condensate is now found for the total impurity 



action Sir 



imp 



imp 



r 

*cav+W drGf mp (-T)* cav , 
Jo 



(A9) 



(A10) 



By plugging this equation into Eq. IA8I a closed set of 
equations for the true condensate is found, which is tan- 
tamount to B-DMFT. The iteration scheme proceeds in 



which can be inverted (accurate up to this order), 
*cav ~ (I - tGf mp (iu) n = 0)) * imp . 



the usual way. Note that the last step does not rely on 
the scaling in Eq. ([75)1 and Eq.([77)l. but does use the re- 
currence relation of a tree. We note that the derivation 
presented in Ref. 1371 is conceptually equivalent to the one 
presented in Ref. |49|, on which this appendix is based. 



Appendix B: Effective medium approach 

In this section we will expand the full Bose-Hubbard 
action around a single site, identify the low-energy 
modes, and re-exponentiate them. The derivation in this 
section is similar to the derivation found in Ref. l33l . but 
differs in the treatment of the condensate. We provide a 
microscopic RG-like picture for the condensate field and 
hybridization fields in the effective action. The total ac- 
tion of the full Bose-Hubbard model is split into three 
parts, with the first one being the local part defined in 
Eq. ©, 

P ^ u 
Shit = I drb* at (r)(d T - fi)b iat (r) + —nint{nint - 1). 

(Bl) 

The second part AS" describes the coupling with the rest 
of the lattice, 



AS = 



[ AS(t)= ( 
Jo Jo 



(int,cxt) 



(B2) 

and the third part is the remaining, external part 'ext' 
with action 

f U 
Soxt = / dTbl xt (T)(d T -fi)b 0Xt + —n cxt (n C xt~l). (B3) 
Jo * 

We had that S = S'int+Soxt+A^, but would like to derive 
an action S ~ Simp + Scxt by expanding the exponential 
function with respect to the action AS, and perform the 
functional integral over the external degrees of freedom, 
and incorporate this into a new effective internal action 
Si m p and a decoupled, external action S cx t- 

We assume that 'ext' is a thermodynamic system that 
can spontaneously break the symmetry, if needed, and 
develop a condensate. If so, we decompose 



Sb. 



(B4) 



where (6 ex t) = 0cxt represents the condensate wave 
function. It is a classical field (c-numbcr, taken to be 
real) that breaks the symmetry, and its presence forces 
us to consider anomalous averages such as (& e xt&ext), 
which are nonzero. The commutation relations are now 
[<5&ext> <^4xt ] = an d zero otherwise. The impurity 
part 'int' is by itself small and cannot spontaneously 
break the symmetry and develop a non-zero expectation 
value for its local operators. 

In our perturbative expansion in AS an atom from the 
'int' part will hop to the 'ext' part and acquire an ex- 
pectation value through the full correlators in 'ext'. A 
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non-zero expectation value corresponds to the conden- 
sate mode, and is not possible for the bosons outside the 
condensate. We do not attempt here to build in all low- 
energy modes via 'ext', only the classical field contribu- 
tions from the condensate will suffice. Let us anticipate 
in advance on such a possibility and allow 



(B5) 



We also have the commutation relations [Sb^.,Sbl^ ] — 

Sj^k and [£&int > *^ext ] = 0- In second order, we then also 
have to allow for (&i n t&int) to become large. The validity 
of Eq. (|B5j) can then be checked a posteriori. 

The coupling between the internal and external degrees 
of freedom, expressed in Eq. (|B2[) . can then be written 



AS 



~ l ( 6 int & e X t + &Lt fo int 

(int,ext) 

p 

dTz4> int (6b cxt + Sbl xt ) + z4> ext (5b h 



Y2 Sb Lt Sb ex* + Sb Lt Sb " 

(int,ext) 



(B6) 



where we have omitted terms in ^int^ext since the conden- 
sate is treated as a c-number and such terms involve an 
arbitrary shift of the condensate energy. There are four 
terms in AS, two of which have the condensate directly 
in them, and two that don't. Since we are interested in 
adding the physics of the WIBG to the action, we have to 
treat the condensate as a large contribution compared to 
the small contributions coming from the depleted atoms, 
while the third term will be contributing to one (and 
higher) loop corrections. We can now immediately write 
down the contribution of the parts with the condensate 
lines to the action, 

Si = zt<p cxt [ Sb(r) + 5b\r) = zt<f> ext f b(r) + b\r) 



(B7) 

and we keep only the following terms in AS (the other 
one factorizes) , which we treat as a small term, 

rP 

A S = -t Sb\ nt Sb cxt + 5bt xt 5b int . (B8) 

^ (int,ext) 

The advantage of our approach is that we have first added 
the condensate to the impurity action, before looking at 
depletion (fluctuation) terms. In this picture, the sta- 
tionarity and time-independence of the condensate is im- 
mediately built in. We are now in a position to expand 
in AS. This results in an infinite series, 

Z ~ Z cxt / D[6* nt ,6 in t]e- 5 '" t l h *»- h '» t ^ Sl [ b *nf b "' t ]C 



rP 

C = 1 - / dr(A5(r)>«t + 
Jo 



i r r 

r / dn dr 2 (A5(r 1 )A5(r 2 )) cxt + . . (B9) 
1 Jo Jo 



We start with the evaluation of the first order term, 

rfi rP 

/ (A5)ext = -W olr V (<5&[ nt £& ext )ext + h.c, 
Jo Jo ,. . ,\ 

(int,ext) 

(BIO) 

which is zero, as could have been expected from symme- 
try. 

The second order term is the lowest order term that 
survives, 

i-P r p 



If dn f dr 2 (A5(T 1 )A5(r 2 )) cxt = 
1 Jo Jo 

%- [ dn [ dr 2 V [S 1 + S 2 + S 3 + S% (Bll) 
2 Jo Jo 



j,fc£cxt 



with S*,j = l,...,4, 



S 1 

s 2 
s 3 
s 4 



= (8b int ( Tl )Sb int (T 2 )Sb^ (T^Sblif ( T2 )) ext , 

= (Sb- mt (n)6bl t ir 2 )5b]^ (n)6bW (r 2 )) cxt , 

= (^t(n)^int(T 2 )^(n)^eit ) (^))ex t , 

= (^Lt(n)^ n t(^)^ ) t(r 1 )<56S(r 2 )) oxt . (B12) 

The notation (j) and (fc) denotes different sites that are 
coupled to the impurity through the hopping term in 
the Bosc-Hubbard Hamiltonian. The anomalous terms 
survive in the presence of a condensate, which tells us 
that there are four terms in second order. Let us analyze 
S 1 in more detail, 

S 1 = ( 56 int (r 1 ) ( 56 int (r 2 )(<5^ x J t ) (r 1 )^ ) (r2))c X t. (B13) 

This describes the anomalous process of a depleted bo- 
son hopping from the impurity site to the environment, 
propagating there according to the full (but unknown to 
us) anomalous two-particle Green's function, and finally 
being annihilated on the impurity site. We will now per- 
form a cumulant re-exponentiation of this term in the 
presence of a condensate, for which only the connected 
diagrams can be kept, 



- /(f dn dT 2 56 lnt (T 1 )K:(T 1 -r 2 )(56 Jnt (T 2 ) 



(B14) 



where K will be a function which we treat as unknown 
and which originates from the connected diagrams only. 
The DMFT selfconsistency iteration scheme will deter- 
mine K selfconsistently. In a numerical algorithm it is 
practical to work with full operators b instead of (56, so 
we will write 6j nt (r) — 4>int instead of 5b lnt . The other 
terms S^,j = 2,3,4 are treated in the same way as S 1 , 
but we will not write this out explicitly. 

After convergence, we wish that our site is in equi- 
librium with the surroundings. We thus impose a first 
selfconsistency relation 



(B15) 
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The meaning of this relation is that in every iteration we We collect now all terms needed for our effective impurity 
have to compute (b) and equate the condensate with this action, 
value. It can always be taken real and time-independent. 

I 



Si 



imp 



/ dTb* nt (T)(d T - /x)6i n t(r) + — n int (n int - 1) - zt*Lt / ^b int (r) 
Jo 1 Jo 

+ \ [ dr £ dr'(bUr) - * int (r))A(r - r')(b int (r') - * int (r')), 



r 



(B16) 



where the elements of the hybridization matrix are (drop- 
ping the subscript 'c' for connected). The remaining hy- 
bridization term can now be treated similarly as the hy- 
bridization term in fermionic DMFT, and has hence a 
similar selfconsistency relation^ 3 , (see Sec. IIV|) . 



Appendix C: High frequency expansions 

Within our selfconsistency loop we need to perform 
Fourier transforms between the imaginary time r domain 
and Matsubara frequencies u n . Even though we use a 
continuous time method we measure the Green's func- 
tion G(t) on a grid with a fixed number of time slices. 
We can therefore only obtain accurate data for a finite 
number of Matsubara frequencies. In order to perform 
accurate Fourier transforms we supplement the Green's 
function with the analytically known high-frequency be- 
havior, given by 



G{iu) n ) = - 



91 



92 



iuj n (iuJ n ) 2 («W n ) 3 



+ AG(iu n ) (CI) 



where AG{iw n ) goes to zero as l/(iw n ) for iuj n — ► 
oo. For the static component (n = 0) we let 
AG(iuj n ) — G(iu! n ). Computing these coefficients ana- 
lytically improves the accuracy of our Fourier transfor- 
mation significantly^ By using 
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1 ^ e - itd " T 
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we can write the inverse Fourier transformation as 
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and the Fourier transformation as 
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+ 



(iw n ) 2 (iu n ) 
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AG(r)e l 



(C5) 

with AG(r) = G(r) - g{r). 

The coefficients for the Green's function are obtained 
by evaluating commutators: 

, 9l = ([6,6+]) = 1, 

,g 2 = ([[H,b],tf])=n-2U(n), 

.93 = ([[#,[#, fo]],&t]) 

= (e 2 ) +fi 2 + 3U 2 (n 2 ) - (n)(4 M £/ + £/ 2 ), 

,9i = ([b,b]) = 0, 

.9 2 = {[[H,b],b]) = U{bb), 

.93 = ([[H,[H,b)},b}) = 0, (C6) 

where gi are the coefficients for the diagonal and gi for 
the off-diagonal Green's function. All odd coefficients of 
the anomalous functions are zero, as these functions are 
purely real since all anomalous functions are symmetric 
in imaginary time. 

Similar coefficients can be obtained for the hybridiza- 
tion function A(iui n ). Here we find 



h = 
h = 
h = 



-<e 2 )52, 

( £ 4 ) + ( e 2 )(.93-2( e 2 )), 



(C7) 



for the diagonal component F(iuj n ) and 

0. 



ki 
k 2 
k 3 



-0.5(e 2 )g 2 



0. 



(C8) 



(C3) for the off-diagonal component K(iui n ). 
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